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ABSTRACT 



A number of experiments are set to measure the 21-cm signal of neutral hydrogen from 
Q the Epoch of Reionization (EoR). The common denominator of these experiments are the large 

\-[ data sets produced, contaminated by various instrumental effects, ionospheric distortions, RFT 

and strong Galactic and extragalactic foregrounds. In this paper, the first in a series, we present 
^ the Data Model that will be the basis of the signal analysis for the LOFAR (Low Frequency 

Array) EoR Key Science Project (LOFAR EoR KSP). Using this data model we simulate real- 
istic visibility data sets over a wide frequency band, taking properly into account all currently 
^ known instrumental corruptions (e.g. direction-dependent gains, complex gains, polarization 

^\ effects, noise, etc). We then apply primary calibration errors to the data in a statistical sense, 

lO assuming that the calibration errors are random Gaussian variates at a level consistent with 

our current knowledge based on observations with the LOFAR Core Station 1 . Our aim is to 
CO demonstrate how the systematics of an interferometric measurement affect the quality of the 

^-H calibrated data, how errors correlate and propagate, and in the long run how this can lead to 

new calibration strategies. We present results of these simulations and the inversion process 
On and extraction procedure. We also discuss some general properties of the coherency matrix 

and Jones formalism that might prove useful in solving the calibration problem of aperture 
^ synthesis arrays. We conclude that even in the presence of reahstic noise and instrumental 

• errors, the statistical signature of the EoR signal can be detected by LOFAR with relatively 

small errors. A detailed study of the statistical properties of our data model and more complex 
H instrumental models will be considered in the future. 
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1 INTRODUCTION 

Recent years have seen a marked increase in the study, both the- 
oretical and observational, of the epoch in the history of our Uni- 
verse after the cosmological recombination era: from the so called 
'Dark Ages' to the Epoch of Reionization (EoR) (Hogan & Rees 
1979; Scott & Rees 1990; Madau, Meiksin, &. Rees 1997). A cold 
and dark Universe, after the recombination era, was illuminated by 
sources of radiation, be it stars, quasars or dark matter annihilation. 
These 'first objects' ionized and heated their surrounding inter- 
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galactic medium (IGM), carving out 'bubbles' in the otherwise neu- 
tral hydrogen-filled Universe. These bubbles grew rapidly, both in 
size and number, and caused a phase transition in the hydrogen- 
ionized fraction of our Universe at redshifts ()<z<2Q (Sunyaev 
& Zeldovich 1975). Although the EoR spanned a relatively small 
fraction, in time, of the Universe's age, its impact on subsequent 
structure formation (at least baryonic) is crucial. Hence, studying 
the EoR directly influences our understanding of issues in con- 
temporary astrophysical research such as metal-poor stars, early 
galaxy formation, quasars and cosmology (Nusser 2005; Zaroubi & 
Silk 2005; Kuhlen & Madau 2005; Thomas & Zaroubi 2008; Field 
1958, 1959; Scott & Rees 1990; Kumar, Subramanian, & Padman- 
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abhan 1995; Madau, Meiksin, & Rees 1997). For a detailed review 
of the EoR and our current efforts to detect it, we refer the reader 
to Furlanetto, Oh, & Briggs (2006) and the references therein. 

Given the recent progress in developing a concrete theoreti- 
cal framework, and simulations based thereon, the EoR from an 
observational point of view is still very poorly constrained. De- 
spite a wealth of observational cosmological data made available 
during the past years (e.g. Spergel et al. 2007; Page et al. 2007; 
Becker et al. 2001; Fan et al. 2001; Pentericci et al. 2002; White 
et al. 2003; Fan et al. 2006), data directly probing the EoR have 
eluded scientists and the ones that constrain the EoR are indi- 
rect and very model-dependent (Barkana & Loeb 2001; Loeb & 
Barkana 2001; Ciardi, Ferrara, & White 2003; Ciardi, Stoehr, & 
White 2003; Bromm & Larson 2004; Ihev et al. 2007; Zaroubi et 
al. 2007; Thomas & Zaroubi 2008). Currently there are two main 
observational constraints on the EoR: first, the sudden jump in the 
Lyman-a optical depth in the Guim-Peterson troughs (Gunn & Pe- 
terson 1965), observed in the quasar spectra of the Sloan Digital 
Sky Survey (SDSS) (Becker et al. 2001; Fan et al. 2001; Pentericci 
et al. 2002; White et al.q 2003; Fan et al. 2006) which provides 
a Umit on when reionization was completed. Current consensus is 
that reionization ended around a redshift of six. Second, the five- 
year WMAP data on the temperature and polarization anisotropies 
of the cosmic microwave background (CMB) (Spergel et al. 2007; 
Page et al. 2007) which gives an integral constraint on the Thomson 
optical depth for scattering experienced by the CMB photons since 
the EoR. A maximum likelihood analysis performed by Spergel 
et al. (2007) estimates the peak of reionization to have occured at 
11.3 when the cosmic age was 365 Myr. Thus, we see that current 
astronomical data is only able to provide us with crude boundaries 
within which reionization occured. In order to properly characterize 
the onset, evolution and completion of the EoR and derive results 
on its impact on subsequent evolution of structures in the Universe, 
we need more direct measurements from the EoR. 

Observations of the hydrogen 21-cm h5'perfine "spin-flip" 
transition, using radio interferometry, provide just such a direct 
probe of the dark ages and the EoR over a wide spatial and redshift 
range. It is worth mentioning that the "spatial range" here implies 
the two dimensions on the sky, which is a function of the baseline 
lengths of an interferometer, and the third dimension along the red- 
shift direction, which depends on the frequency resolution of the 
observation. The 21-cm emission line from the EoR is redshifted 
by 1 -h because of the expansion of the Universe, to wavelengths 
in the meter waveband. For example, at a redshift of « = 9 the 
21-cm line is redshifted to 2.1 metres, which corresponds to a fre- 
quency of about 140 MHz. Computer simulations suggest that we 
may expect a complex, evolving patch-work of neutral (Hi) and 
ionized hydrogen (Hn) regions. If we manage to successfully im- 
age the Universe at these high redshifts (6 < 2; < 12) we expect 
to find Hii regions, created by ionizing radiation from first objects, 
to appear as "holes" in an otherwise neutral hydrogen-filled Uni- 
verse; the so-called Swiss-cheese model. Current constraints and 
simulations converge on reionization happening for a large part in 
the redshift range z w 11.4 (~115 MHz) to z w 6 (~203 MHz), 
which is the range probed by LOFAR^ , a radio interferometer cur- 
rently being built near the village of Exloo in the Netherlands. The 
21-cm radiation can not only trace the matter power spectrum in the 
period after recombination, but also can constrain reionization sce- 
narios (Thomas & Zaroubi 2008; Barkana & Loeb 2001). Note that 
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because strongly radiating sources create bubbles of ionized gas in- 
side the neutral IGM, one should observe fluctuations in the 21-cm 
emission due to reionization that deviate from those of neutral gas 
tracing the dark-matter distribution, even deep into the highly linear 
regime. 

Developments in radio-wave sensor technologies in recent 
years have enabled us to conceive of and design extremely large, 
high sensitivity and high resolution radio interferometers, a devel- 
opment which is essential to conduct a successful 21-cm experi- 
ment to image the EoR. A series of radio telescopes are being built 
similar to LOFAR, such as MWA^ PAPER ^ 21CMA^ and furtiier 
in the future the SKA®, all with one of their primary goals being the 
detection of the redshifted 21-cm signal from the EoR. The GMRT^ 
has a programme akeady under way to detect the EoR or at least to 
constrain the foregrounds that may hamper the experiment (Pen et 
al. 2008). 

Calculations predict the cosmological 21-cm signal from the 
EoR to be extremely faint. Apart from the intrinsic low strength 
of the 21-cm signal, the experiment is plagued by a myriad of 
signal contaminants like man-made and natural (e.g. lightnings) 
interference, ionospheric distortions. Galactic free-free and syn- 
chrotron radiation, clusters and radio galaxies along the path of 
the signal. Thus long integration times, exquisite calibration and 
well-designed RFI mitigation techniques are needed in order to 
ensure the detection of the underlying signal. It is also impera- 
tive to properly model all these effects beforehand, in order to de- 
velop sophisticated schemes that will be needed to clean the data 
cubes from these contaminants (Shaver et al. 1999; Di Matteo, Cia- 
rdi, & Miniati 2004; Di Matteo et al. 2002; Oh & Mack 2003; 
Cooray & Furlanetto 2004; Zaldarriaga, Furlanetto, & Hemquist 
2004; Gleser, Nusser, & Benson 2007). Due to the low signal-to- 
noise ratio per resolution element (of the order of 0.2 or for LO- 
FAR and and even less for e.g. MWA), the initial aim of all current 
experiments is to obtain a statistical detection of the signal. By sta- 
tistical, we mean a global change in a property of the signal, for 
example the variance, as a function of frequency and angular scale. 
Note that this task involves distinguishing these statistical proper- 
ties from those of the calibration residuals and the thermal noise. 

In addition to the above-mentioned astrophysical and terres- 
trial sources of contamination, one also has to face issues arising 
in standard synthesis imaging. For that it is crucial to describe all 
physical effects on the signal that determine the values of the mea- 
sured visibiUties. The study of polarized radiation falls within the 
regime of optics: Hamaker & Bregman (1996) provided such a 
unifying model for the Jones and Mueller calculi in optics (Bom 
& Wolf 1999) and the techniques of radio interferometry based 
on multiplying correlators. Because low frequency phased-array 
dipole antennas are inherently polarized, one has to consider po- 
larimetry from the beginning. The Measurement Equation (ME) 
of Hamaker & Bregman (1996) is therefore a natural way to de- 
scribe LOFAR. Their treatment (Hamaker & Bregman 1996; Sault, 
Hamaker, & Bregman 1996; Hamaker 2000c, 2006) forms the ba- 
sis of our data model description and we will present it in a wider 
setting, giving the connection to physics. The Hamaker-Bregman- 
Sault measurement equation acts on the astronomical signals as 

Murchison Widefield Array: http://www.haystack.mit.edu/ast/arrays/mwa 
Precision Array to Probe Epoch of Reionization: 
http://astro.berkeley.edu/ dbacker/eor/ 

21-cm Array: http://web.phys.cmu.edu/ past/ 
^ Square Kilometer Array: http://www.skatelescope.org 
^ Giant Meterwave Radio Telescope: http://www.gmrt.ncra.tifr.res.in 
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Figure 1. Left: The LOFAR compact-core layout situated near Exloo in the Netherlands. The open circle-size corresponds to the High Band Antenna (HBA) 
station size of approximately 31 metres in diameter. The 6 inner stations ("six-pack") are shown in the the inset figure. Middle: The snapshot uv-coverage of 
the LOFAR compact core at the zenith and at a frequency 150 MHZ. Riglit: The layout of a LOFAR HBA station. Each square represents a tile that consists 
of four by four, orthogonal bowtie dipoles, as shown in the blown-up inset. Each dipole pair can be rotated individually during the installation of the tile. 



a "black-box": the interferometer converts the input signal Stokes 
vectors to the final output at the correlator. This is done via a se- 
quence of linear transformations and thus enables us to systemat- 
ically model the series of effects that modify the signal while it 
propagates through the ionosphere and the receivers. 

Calibration (Morales & Matejek 2008) of the observed visibil- 
ity data set is generally aimed at determining instrumental parame- 
ters of the antennas at a level sufficient to detect signals at several 
times the noise level. While this is the traditional approach, a much 
more thorough understanding of the instrumental response is re- 
quired in currently designed experiments such as the LOFAR EoR 
Key-Science Project, where unprecedented high dynamic range and 
sensitivity have to be achieved, and the signal is far below the noise 
level. For example, the CLEAN algorithm and variations on it have 
been used extensively, as an integral part of the SELFCAL pro- 
cess (Pearson & Readhead 1984). While computationally efficient, 
it does not provide a statistically optimal solution (Schwarz 1978; 
Starck, Pantin, & Murtagh 2002). In this work we shall therefore 
consider a maximum-likelihood solution to the measurement equa- 
tion inversion problem (Boonstra 2005; Leshem, van der Veen, & 
Boonstra 2000; Leshem & van der Veen 2000; Wijnholds, Breg- 
man, & Boonstra 2004), which takes into account realistic zero- 
mean calibration residuals and noise. We are also examining alter- 
native solutions (Yatawatta et al. 2008), however. 

After giving a short description of the LOFAR array, with em- 
phasis on the high band (HBA) aspects of the design in Section 2, 
we review the basic relation between the observed visibilities and 
the sky intensity in section 3. Wc emphasize the polarized, matrix 
formulation of the measurement equation and the mathematical as- 
pects of coherency matrix (Hamaker & Bregman 1996). In Section 
4 we briefly discuss the relevant Measurement Equation (ME) pa- 
rameters for LOFAR. Using this measurement equation as our data 
model, in Section 5, we produce a number of simulations for dif- 
ferent instrumental parameters. This is the forward use of the data 
model. We also try to invert the data model, given the instrumen- 
tal parameters and their error distributions, in order to recover the 
original data. Our goal is to test the calibration and inversion re- 
quirements using reaUstically generated data cubes. In Section 6 
we discuss the results and give our conclusions. 



2 DESCRIPTION OF THE LOW-FREQUENCY ARRAY 

The immediate science goals of the LOFAR EoR-KSP, that drive 
some of the considerations about the design of LOFAR (Bregman 
2002), are: (1) extract the 21-cm neutral hydrogen signal averaged 
along lines of sight, i.e. the 'global signal' (e.g. Shaver et at. 1999; 
Jelic et al. 2008), (2) determine the spatial-frequency power spec- 
trum of the brightness temperature fluctuations on angular scales 
of about 1 arcminute to 1 degree and frequency scales between 
0.1 and 10 MHz in the redshift range of ~6-ll and (3) search for 
Stromgren ionization bubbles around bright sources and the 21-cm 
absorption-line forest (de Bruyn et. al 2007). In order to achieve 
these science goals, LOFAR requires a good uv-coverage, a good 
frequency coverage and a large collecting area. 

Below, we give a short summary of the aspects of LOFAR 
which are relevant to the EoR experiment and our data model. For 
more details we refer to the project paper by de Bruyn et al. (in 
preparation). 

2.1 Station configuration and uv-coverage 

In its current layout, the LOFAR telescope (de Bruyn et. al 2007; 
Falcke et al. 2007) will consist of up to 48 stations of which ap- 
proximately 24 will be located in the core region (Figure 1), near 
the village of Exloo in the Netherlands. The core marks an area 
of 1.7 by 2.3 kilometres. Each High Band Antenna station (HBA 
station; 110-240 MHz; see next section) in the core is further split 
into two "half- stations" of half the collecting area (~31 metre di- 
ameter), separated by ~ 130 metres. This split further improves the 
uv-coverage, though at the cost of quadrupling the BlueGene/P cor- 
relator demands (i.e. there are four times more base-lines). The 
central region of the core consists of six closely-packed stations, 
"the six-pack", to ensure improved coverage of the shortest base- 
lines necessary to map out the largest scales on the sky, such as 
the Milky Way. The station-layout yields a snapshot uv-coverage 
at the zenith as shown in Figure 1 . The uv-coverages for a typical 
synthesis time of 4 hrs are shown in Figure 2. 

A good uv-coverage is crucial for several reasons. First to 
improve sampling of the power spectrum of the EoR signal (San- 
tos, Cooray. & Knox 2005; Hobson & Maisinger 2002; Bowman, 
Morales, & Hewitt 2008, 2006, 2005; Morales, Bowman, & He- 
witt 2005). Second, to obtain precise Local (Nijboer, Noordam, & 
Yatawatta 2006) and Global (Smknov & Noordam 2004) Sky mod- 
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Figure 2. The first column of figures shows the uv-coverage of the LOFAR core, that will be used for the EOR experiment, for diffrent values of the declination 
5. The size of each point con'esponds to the station diameter. The uv coverage is calculated for 4 hours of synthesis with 10s averaging at 150 MHz. The 
second column shows the coiTesponding uv point density. The uv plane is gridded with a cell size of 8.5 wavelengths squared. In the last column a horizontal 
cut of the "dirty" beam is shown. 



els (LSM/GSM; i.e. catalogues of the brightest, mostly compact, 
sources in and outside of the beam, i.e. local versus global). A fur- 
ther complication is the extraction of the Galactic and extragalactic 
foregrounds. This is a vital step in the recovery of the signal and 
requires good sampling of the uv-plane at all frequencies (Bow- 
man, Morales, & Hewitt 2008, 2006; Morales, Bowman, & Hewitt 
2005). 



2.2 The High-Band Antennas 

LOFAR will have two sets of dipoles, the Low Band Antennas 
(LBA) and the High Band Antennas (HBA). For the EoR experi- 
ment we are mostly interested in the HBA dipoles which cover the 
110 to 220 MHz frequency range. Each dipole is a crossed dipole 
which enables X and Y polarization observations. Inside the sta- 
tion the dipoles are arranged in tiles of 4 by 4 dipoles, with 24 tiles 
per station inside the core (Figure 1). Radio waves are sampled 
with a 16-bit analog-to-digital converter to be able to cope with ex- 



pected interference levels operating at either 160 or 200 MHz in the 
first, second or third Nyquist zone (i.e., 0-100, 100-200, or 200- 
300 MHz band respectively for 200 MHz sampling). The data from 
the receptors are filtered in 512, 195 kHz sub-bands (156 kHz sub- 
bands for 160 MHz sampling) of which a total of 32 MHz band- 
width (164 channels) can be used at one time. Sub-bands from all 
antennas are combined at the station level in a digital beamformer 
allowing multiple (4-6) independently steerable beams, which are 
sent to the central processor via a glass-fibre link that handles 0.7 
Tbit/s data. The beams from all stations are further filtered into 1 
kHz channels, cross-correlated and integrated. The integrated visi- 
bilities are then calibrated on 1 second intervals, to correct for the 
effects of the ionosphere, and subsequently images are produced. 
Channels with disturbing radio frequency interference (RFl) are ex- 
cised (Leshem, van der Veen, & Boonstra 2000; Veen, Leshem, & 
Boonstra 2004b; Wijnholds, Bregman, & Boonstra 2004; Fridman 
& Baan 2001). For the correlation we use three racks of an IBM 
Blue Gene/P machine in Groningen with a total of 12288 process- 
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ing cores. LOFAR is a new concept in array design, a broad-band 
aperture array with digital beamforming. This malces LOFAR es- 
sentially qualify as a pathfinder for the Square Kilometre Array 
(Falcke et al. 2007). 



3 MATHEMATICAL FRAMEWORK OF THE LOFAR 
DATA MODEL 

The most important part of any physical measurement is to find a 
correspondence between the physical quantities and the measured 
quantities. In radio interferometry this is achieved through the so 
called measurement equations (Hamaker & Bregman 1996; Boon- 
stra 2005). The measurement equation (ME) describes the relation- 
ship between the visibilities (correlations between the electric fields 
from different antennas) and the brightness distribution of the sky. 
We will begin by discussing briefly the different types of measure- 
ment equations and the implications for LOFAR (Smimov & No- 
ordam 2006; Noordam 2004, 2000). The data model presented in 
this paper can be easily applied to other telescopes that operate at 
low frequencies such as the MWA and eventually the SKA. 

Astronomical radio signals appear as spatially wide-band ran- 
dom noise with superimposed features, such as polarization, emis- 
sion and absorption lines. The physical quantity that underlies this 
kind of measurement is the electric field, but for convenience as- 
tronomers try to recover the intensity in the direction of the unit 
pointing vector s, //(s) = (|S/(s)|)^. The measured correlation 
of the electric fields between two sensors i and j is called the com- 
plex visibility. For Earth-rotation synthesis we assume that the tele- 
scopes have a small field of view (FOV) and that they track a posi- 
tion on the sky. To achieve that, a slowly time-varying phase delay 
has to be introduced at the receiver to compensate for the geometri- 
cal delays. The result is that the reference location appears to be at 
the zenith, or their chosen point in the sky (phase reference center). 
For a planar array, the receiver baselines can be parametrized as 

Ti-Tj = X[U,V,0], \=^. (1) 

where ri are the station position vectors. This system is wavelength 
dependent. The (scalar) measurement equation in (w, v) coordi- 
nates becomes 



Vf (u,v) = jj Vf{Lm)If{l,m)c-''-'''+'"'^dldm 



(2) 



where m) is the complex primary beam or antenna response 
pattern and //(/, m) is the sky brightness distribution. This equa- 
tion (van Cittert-Zernike theorem) is in the form of a 2-D Fourier 
transform, which is an approximation for a flat sky (Thompson, 
Moran, & Swenson 2001; Carozzi & Woan 2008). The visibilities 
are sampled for all different sensor pairs i and j, but also for dif- 
ferent sensor locations projected on the sky, since the Earth rotates. 
Hence Earth-rotation synthesis traces uv tracks for each baseline 
(Thompson, Moran, & Swenson 2001). 

3.1 The Scalar Measurement Equation 

The most widely used ME is still the scalar formulation of the ME. 
We begin the discussion with the scalar ME and later we will also 
discuss the polarized version thereof. Since all processing is done 
with digital computers this equation must be transformed into a 
more convenient discretized form. The main output of the LOFAR 
correlator is a set of correlation matrices (Boonstra 2005; Falcke 
et al. 2007), R/(tfc), for a set of narrow-band frequency channels 



(adding up to 32 MHz bandwidth) and for a set of short-time in- 
tegrations (adding up to >300 hours of integration). The connec- 
tion between the correlation matrices R/ (tk) and the visibiUties 
Vf {u, v) is that each entry Rij (tk) of R/ (tk) is a sample of the 
visibility function for a specific coordinate (u, v) corresponding to 
the baseline vector Vij = — between telescopes i and j at time 
tk (Boonstra 2005): 

Vf [uik - Ujk,Vik - Vjk) = Rij (tk) . 

The noiseless scalar measurement equation (not accounting for in- 
strumental and other distorting effects) for one short-time integra- 
tion and narrow-band frequency channel, assuming the sky can be 
described by a set of point- sources, can then be written in terms 
of the correlation matrices^ as (Boonstra 2005; Leshem & van der 
Veen 2000) 



R 



(3) 



where 



Afc,/ = [afc,/(si), . . . ,afc,/(sd)] 



afc,/(s<) 



Bf = 



-z(uifc(+-uifem) 



-<"«tel'='+''«tel'='") 



(4) 



Bf (si) 



Bfis,) ) 



and where k is the time-ordered visibility number, s is the source 
position vector on the celestial sphere and / the observing fre- 
quency channel number. 

The vector functions afc j are called array response vectors in 
array signal processing and they are frequency dependent, but also 
time dependent in this case due to the rotation of the Earth (Leshem 
& van der Veen 2000; van Trees 2002). They describe the response 
of an interferometer to a source at direction s = {l,m) (see Figure 
3). 

The above formalism is trivial as long as the positions of the 
telescopes are well known. In reality though, the response of the 
array is not perfect: telescopes are not omni-directional antennas, 
but each one has its own properties (i.e. complex beam-shape and 
gain etc.). In this case the array response vectors must be redefined 
as 

-i(uifc-Si) 



afc,/(si) = 







(5) 



where indicates a Hadamard (element-wise) product. The source 
structure can also vary with frequency. Finally, most of the received 

signal consists of additive noise. When the noise has zero mean and 
is independent among the antennas (spatially white), then 

Rfcj = AkjBAlj + CT/I 

Noise can be assumed to be Gaussian in radio-interferometers, like 
LOFAR. We will assume this in the remainder of this paper, and 
might address non-Gaussian (Thompson, Moran, & Swenson 2001) 
time-varying noise in future publications. Actually, system noise is 



^ The symbol "f" stands for the Hermitian conjugation operator 
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slightly different at each receiver. It is reasonable to assume that 
noise is spatially white: the noise covariance matrix is diagonal. 

3.2 The Polarization Measurement Equation 

The equations above describe the relation of the visibilities to the 
total intensity of the source, i.e. the classical Stokes /. To take po- 
larization into account, several modifications are required. The use 
of the polarization information is of great importance for several 
reasons. First, this information provides insight in the physical pro- 
cesses that might exist in the astronomical object of interest. Sec- 
ond, modem telescopes like LOFAR are also inherently polarized. 
Third, using more information enhances the result of the data re- 
duction process. The study of polarized light is therefore becoming 
an increasingly important issue in astrophysics (Tinbergen 1996). 

Many matricial models have been developed to study the po- 
larization properties of light. A proper description of the polar- 
ization properties of light relies on the concept of the coherency 
matrix (Born & Wolf 1999; Hamaker & Bregman 1996). This 
mathematical formalism holds for every band of the electromag- 
netic spectrum. A usual assumption is monochromatic light, but 
polychromatic light behaves as monochromatic for time intervals 
longer than the natural period and shorter than the coherence time 
(Gil 2007; Barakat 1963). Our mathematical model will be based 
on those introduced in radio-astronomy by Hamaker & Bregman 
(1996). 



3.2. 1 The Electric-Field Vector and Coherency Matrix 

The effects of linear passive media on the propagated photons can 
be represented by linear transformations of the electric field vari- 
ables. The nature of those effects, the spectral profile of the light 
and the chromatic and polarizing properties of the medium through 
which light passes, all affect the degree of mutual coherence. In 
general coherent interactions can be represented by the Jones cal- 
culus (Jones 1941, 1942, 1948), while incoherent interactions of 
polychromatic light require the Mueller calculus (Barakat 1963), 
since the loss of coherence needs more parameters to be described. 
The two components of the electric field (e.g. those received at two 
dipoles) can be arranged as the components a 2 x 1 complex vector: 



e(t) = 



Ey (t) e"*(') 



where 5 (t) is the relative phase. This vector includes all informa- 
tion about the temporal evolution of the electric field. When the 
parameters have no time dependence this is called the Jones vector. 
Moreover, the coherency (or polarization or density) matrix of a 
light beam contains all the information about its polarization state. 
This Hermitian 2x2 matrix is defined as 



{e2{t)el{t)) {e2(t)e^(t)) 



(7) 



This is the coherency matrix of the perpendicular dipoles of a sin- 
gle LOFAR HBA antenna, (g) stands for the Kronecker product and 
the brackets indicate averaging over time (Boonstra 2005). The co- 
herency matrix is a correlation matrix whose elements are the sec- 
ond moments of the signal. Using the ergodic hypothesis the brack- 
ets can be considered as ensemble averaging. Due to its statistical 
nature its eigenvalues ought to be non-negative. The normalized 



version of this matrix C = 



tr{C) 



contains information about the 



population and coherencies of the polarization states (Fano 1957). 
This object is the equivalent of the single brightness point in the 
scalar version of the theory. 

3.2.2 The Stokes Parameters 

Above we discussed the statistical interpretation of the coherency 
matrix. Now we can also introduce a geometrical description: the 
measurable quantities Stokes I, Q, U and V arise as the coeffi- 
cients of the projection of the coherency matrix onto a set of Hermi- 
tian trace-orthogonal matrices, the generators of the unitary SU(2) 
group plus the identity matrix (see Appendix Al). Parameters with 
direct physical meaning can be derived from the corresponding 
measurable quantities. The Stokes parameters are usually arranged 
as a 4 X 1 vector, 

f ' \ 

Q 

u 

V V J 

An alternative notation, that will be derived in Appendix A is the 
2x2 Stokes matrix: 



(6) 



s = i 

2 



I + Q 
U + iV 



U-iV 
I-Q 



(8) 
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which relates the measured coherency matrix quantities to the 
Stokes parameters (Bom & Wolf 1999). 



3.2.3 The Jones Formalism 

An adequate method to describe a non-depolarizing system is the 
Jones formalism. It represents the effects on the polarization prop- 
erties of an EM wave after the interaction with such a system. For 

passive, pure systems, the electric field components of the light in- 
teracting with them is given by the corresponding Jones matrix J, 

e' — Je. 

As both the initial and final fields can fluctuate, it is useful to de- 
scribe the properties of partially polarized light with the coherency 
matrix. Thus, 



C = ^e' ® e't) = ^(Je) ® (Je)+J> 



(9) 



As we are dealing with interferometry, the two J matrices can 
come from two different telescopes. The effects on the elec- 
tric field vectors in the coherency matrix C' can be written 
as an operation of these Jones matrices on the original unaf- 
fected coherency matrix C. As we already mentioned, the co- 
herency matrix can also be written as a four vector with c = 
((eie^) , (6162) , (e2et) , (6262)). This vector is related to the 
Stokes vector via (Parke 1948) 



s = Lc with L : 




(10) 



The matrix has the following property = V2-'-'^- Using the 
properties of the Kronecker product we then find, in terms of Stokes 
parameters, that 



s' =L(Je® (Je) 



: NS 

N = L(J(8)J)L"' 
with 



(11) 



Nfc; = ^tr ( (JkJcTiJ 



where ctj are the Pauli matrices (Appendix Al). A Jones matrix can 
represent a physically realizable state as long as the transmittance 
condition (gain or intensity transmittance) holds; that is, the ratio of 
the initial and final intensities must be ^ p ^ 1. The reciprocity 
condition describes the effect when the output signal follows the 
path in the inverse order. For every proper Jones matrix e' = J^e. 
This result does not hold when magneto-optic effects are present. 
In this case the Mueller-Jones matrices have to be used. If a Jones 
matrix represents a physically realizable state the reciprocal matrix 
also represents a physical effect. 



4 THE LOFAR MEASUREMENT EQUATION 

LOFAR stations consist of tiles of 4 x 4 crossed dipoles which allow 
for full Stokes measurements. If we have N^ei stations each with 
two polarization degrees of freedom, then the 2Ntei electric fields 



can be stacked into a single vector and the correlation matrix can 
be generalized to the following form 



{Ey{r{)El{vi)) 



E.{vi)E;{r,)) 
Ey in) E; {rj)) 



(12) 



for a linear polarization basis. Note that each element of V, Vij is 
not Hermitian for i ^ j, but Vij = V^^, so that V remains Hermi- 
tian (Hamaker & Bregman 1996). Here we have denoted explicitly 
the correlation from the X and Y oriented dipoles. 

Let Ail be the position dependent polarization multiplication 
matrix. This is the array response matrix (as it is termed in the 
language of communication theory) or in the Hamaker formalism 
the Jones matrix. The array response vector must take into account 
the different physical and instrumental effects that affect the signal 
through its path from the source to the recorder, like ionospheric 
Faraday rotation, parallactic offsets, the geometric delay and in- 
strumental gains (and leakages). These effects are represented by 
the Jones matrices. So we can write the measurement equation in 
the form (Boonstra 2005; Veen, Leshem, & Boonstra 2004a) 

Yiji = AiiB;Aj, -l-Nnoise 

where the different effects A in the array response vector must be 
introduced in the exact order in which they affect the signal. The 
index / is the pixel number and the index i, j represent antenas i 
and j respectively The algebra of complex Jones matrices is obvi- 
ously non-commutative, i.e. the ordering of the matrices matters. 
The physical meaning of this is that the results of the different ef- 
fects on the incoming electromagnetic wave are not Unear. 



4.1 Individual Jones matrices 

In this section we give a brief overview of the instrumental parame- 
ters (Jones matrices) which are used for our data simulations, their 
importance, and which parameters we use and perturb. 

F: lonosplieric Faraday Rotation The ionosphere is birefringent, 
such that one handedness of circular polarization is delayed with 
respect to the other, introducing a dispersive phase shift in radians 

A</> « 2.62 X 10""A^ j B||neds(SI) 

It rotates the linear polarization position angle and is more impor- 
tant at longer wavelengths, at times close to the solar maximum and 
at Sunrise or Sunset, when the ionosphere is most active and vari- 
able (EUasson & Thide 2008; Norin et al. 2008; Eliasson & Thide 
2007; Thide 2007; Spoelstra 1996, 1995; van Velthoven & Spoel- 
stra 1992; Spoelstra & Yang 1990; Spoelstra & Schihzzi 1981). The 
Faraday rotation Jones matrix is a real-valued rotation matrix: 



cos A0 
sin 



— sin A(/i 
cos A0 



To model the Total Electron Content (TEC) (van der Tol & van der 
Veen 2007), we use a Gaussian random field with Kolmogorov- 
like turbulence. We use a cut-off at scales that correspond to the 
maximum and minimum baseUnes of WSRT (Westerbork Synthe- 
sis Radio Telescope), similar to those of LOFAR. Initial model pa- 
rameters have been obtained from WSRT data, taken at similar fre- 
quencies, baselines and time intervals to our planned LOFAR ob- 
servations. Moreover, the WSRT is located approximately 50 km 
from the LOFAR core in Exloo and as such the WSRT provides a 
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Figure 4. The design of a LOFAR HBA dipole antenna element. 




good test-bed for forthcoming LOFAR observations and expected 
ionospheric effects. 

P: Parallactic Angle This matrix describes the orientation of the 
sky in the telescope 's field of view, or similarly the projection of 
the dipoles onto the sky. It has the mathematical structure of a ro- 
tation matrix and rotates the position angle of linearly polarized 
radiation incident on the dipoles. It should in principle be analyti- 
cally and deterministically known, and its variation provides lever- 
age for determining polarization-dependent effects. It can also be 
used at the initial stages to determine dipole orientation errors. 

cosx -sinx 
sin X cos X 

where x is the parallactic angle. We shall include this matrix in the 
antenna voltage pattern matrix. 

E: Antenna Voltage Patterns The LOFAR telescope HBA sta- 
tions consists of crossed bowtie dipole pairs. Like every antenna, 
they have a directionally dependent gain which is important when 
the angular size of the region on the sky is comparable to ~ X/D, 
where D is the station diameter. At low frequencies when the ra- 
dio sky is dominated by point sources, wide-field techniques are 
required as well. To determine the antenna voltage pattern to first 
order, we assume an analytic model of the dipole pairs: the basic 
parameters as shown in Fig. 4 have the following values: L=0.366 
m, h=0.45 m, qi= 50 deg and a2 = 80 deg. 

As the Earth rotates, the dipoles rotate with respect to the sky. 
This causes the polarization coordinate system to rotate (unlike e.g. 
the WSRT which has an equatorial mount) and we must thus correct 
for this effect. For a pair of crossed dipoles along the X and Y axes 
the relevant Jones matrix E can be written as: 



4 / 
4 ) 



(f 



4 / 
3?r "I 
4 ) 



where (p and are the polar and azimuthal angles of the beam pat- 
tern. We assume that the X and Y dipoles have the same beam 
pattern. Because of the spatial distribution of the dipoles within 
a station, each dipole has a different delay for the reception of an 
electromagnetic wave coming from a source. If we assume a nar- 
rowband system, one can correct for this by shifting the phase for 
the signal of each dipole element of the stations (i.e. introducing an 
effective delay). We have an analytical model for the HBA antenna 



element beam and we assume that the dipoles are uniformely dis- 
tributed to form a circular array. Thus, we ignore the tile structure, 
for the purposes of the current paper. If we need to form a station- 
beam in an arbitrary direction inside the dipole/ tile beam, we can 
do it by properly weighting the signals of each element and adding 
delays. Another simpler way to proceed, which is sufficient for the 
purposes of this paper, is to calculate the beam shape of an axially 
symmetric, uniformly distributed array of elements and then mul- 
tiply it with the primary element pattern that we calculated above. 
We also assume that the thin-wire approximation holds. 

D: Polarization Leakage and Instrumental Polarization Radio- 
interferometers usually measure the full set of Stokes parameters 
simultaneously. Measuring the polarization properties of Galactic 
diffuse emission as well as the extragalactic sources is an impor- 
tant aspect of the EoR-KSR The polarization fraction of most as- 
tronomical sources is typically low, of the order of few per cent. 
Measuring it with accuracy is thus challenging, but important in 
order to extract scientific information from it. Cross-polarization or 
polarization leakage is an instrumental contamination that couples 
orthogonal Jones vectors. The dipoles are not ideal, so orthogonal 
polarizations are not perfectly seperated. A slight tilt between the 
dipoles, for example, can change the polarization reception pattern. 
This must be taken into account, especially when we aim to reach a 
high dynamic range in the images. It is supposed to be of the order 
of less than a few per cent and because it is a geometrical factor, 
we expect that it is frequency dependent. The relevant Jones matrix 
can be written in the form of a unitary matrix (Heiles et al. 2001; 
Heiles 2002; Bhatnagar & Nityananda 2001; Reid et al. 2008): 



D 



- iip 



(13) 



where d is a generally small constant of the order of 10~^ This 
is a first-order approximation to behaviour that might prove more 
complex. However, it is very important to know this parameter as 
it can convert unpolarized radiation, such as the EoR signal, into 
polarized radiation. This matrix has almost the structure of a rota- 
tion matrix. The difference is the inclusion of the phase tp. When 
t/j equals zero, polarization leakage converts Q to U or i5 to B- 
modes [rotation across the equator of the Poincare sphere (Heiles 
et al. 2001)]. When the phase term is significant, polarization leak- 
age leads to mixing Q, U and V. Polarization leakage manifests 
itself as closure errors in parallel hand visibilities and the leakage- 
induced closure phase and the Pancharatnam phase of optics (Berry 
1987). This phase arises when the state of polarization of the light 
is transformed following a closed path in the space of states of po- 
larization, which is known as the Poincare sphere (Born & Wolf 
1999). Another type of distortion is the instrumental polarization. 
The dipole and station designs as well as the processing done can 
contribute to such effects. LOFAR uses linearly polarized dipoles. 
The interferometric measurement of the Stokes parameters U and V 
using two different dipoles is made by forming the "cross-handed" 
products of the signals. If for any reason the signals are received 
with different gains, this would lead to a certain fraction of polar- 
ization. This can be modelled as: 



di 






d2e''t' 



(14) 



where di, d2 are the different X and Y gains. We assume that the 
instrumental polarization axis lies along the axis defined by the ma- 
trix basis to write this equation. In such case Q does not leak into 
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U, but still / leaks onto Q. With a carefully chosen basis (see Ap- 
pendix A) / does not leak into U. 

G: Complex baseline-based electronic gain This matrix ac- 
counts for all amplitude-phase- and frequency-independent effects 
of the station electronics. It shows a slow variation of the order of 
1-2% over a day. It has the form: 



G 



gx 

gy 



where gx and gy are complex numbers. It is one of the most com- 
mon calibration parameters and the one commonly solved for in the 
classical self-calibration loop. 

B: Bandpass Compensating for change of gain with frequency is 
called bandpass calibration. This matrix is similar to G and de- 
scribes the frequency-dependence of the antenna electronics. Sta- 
tion digital poly-filters, used to select the frequency passband, are 
not perfectly square, but they are deterministically known and their 
bandpass behavior is expected to be quite stable over large periods 
of time. Its generic form is: 

b.if) 

by if) 

The b are real numbers. We model both B and G parameters as 
slowly time-varying functions with additive white noise as the ini- 
tial comissioning tests for LOFAR suggest. 

K: Fourier kernel This term has traditionally been called the 
Fourier kernel. For larger fields of view and/or longer baselines the 
2-D Fourier transform relationship between the sky and the mea- 
sured visibilities is no longer a good approximation. This term de- 
scribes a phase shift that accounts for the geometrical delay for the 
signals received by telescopes i and j. Let x, y and z be the antenna 
positions in a coordinate system that points to the West, South and 
the Zenith respectively. We can provide the directions on the sky in 
terms of the directional cosines I, m, n. We assume that the sky is a 
unit sphere so that 1^ + + ra^ = 1 (Thompson, Moran, & Swen- 
son 2001). Since the X and Y dipoles are co-located, the relevant 
Jones matrix takes the following simple form: 



K 



— i ^ul-\-v7n-\-iv i^yj 1— i^— m^ — 1^ j 



1 
1 



This matrix should also be well determined, but solving for it can 
help to estimate the antenna positions, the electtonic path lengths, 
the clock errors etc. For this the concept of the array manifold 
is useful (Manikas 2004) It also provides astrometric information 
by solving for Z, m and n, which are the direction cosines of the 
sources on the sky. 

4.2 The Final Form of the Measurement Equation 

All the above effects can be combined in order to form the 
Hamaker-Bregman-Sault measurement equation. The equation 
looks hke: 

V,f^ = f dldm (^Ai (gi AJ^ c (/, m), with 



Ai — KiBiGiDiEiPiT^Fi 



(15) 



where c = vec(C) is the 4x 1 source coherency vector in the 
l,n,m system*. 

If we assume that the sky can be described as a closely-packed 
collection of point-sources, then the measurement equation can be 
rewritten as: 

Fif " = X! Ai (Z, m) C {I, m) A] (Z, m) 

where C is the coherency matrix of a single point-source at {I, rn). 
For a single point-source with the rest of the sky blank this simpli- 
fies further to 

V,f' = Ai {I, m) C (Z, m) A+ {I, m) . 

We must note that this equation is linear over the coherency matrix 
C, but not over the Jones matrices A. We reiterate that the invidual 
Jones matrices that form the matrix products A, are not in general 
commutative. Their order follows the signal path. The form of this 
equation is quite complicated. Every element in the new matrix is a 
non-Unear function of all the parameters that appear in the various 
Jones matrices. 

It has been proposed (Hamaker 2006, 2000a) that we can self- 
calibrate the generic A matrix, apply post-calibration constraints to 
ensure consistency of the astronomical absolute calibrations, and 
recover full polarization measurements of the sky. This is an in- 
teresting idea for low-frequency arrays, where isolated calibrators 
are unavailable (due to the fact that such arrays see the whole sky). 
This just emphasizes the fact that calibration and source structure 
are tied together - one cannot have one without the other. 



4.3 Additive errors 

The ultimate sensitivity of a receiving system is determined prin- 
cipally by the system noise. The discussion of the noise properties 
of a complex receiving system like LOFAR can be lengthy (Lopez 
& Fabregas 2002; de Brujm et. al 2007), so we concentrate for our 
purposes on some basic principles. The theoretical rms (root mean 
square) noise level in terms of flux density on the final image is 
given by (Thompson, Moran, & Swenson 2001) 



SEFD 



Vs X {N -1) X Ai/x ti, 



(16) 



Al= FtTtptEtDtG^BTKT 

J 3 3 3 3 3 3 3 3 



itx.'t 



where tJs is the system efficiency that accounts for electronic, dig- 
ital losses, A'^ is the number of substations, Aj/ is the frequency 
bandwidth, iint is the total integration time and SEFD is the Sys- 
tem Equivalent Flux Density. The system noise we assume to have 
two contributions: the first comes from the sky and is frequency de- 
pendent (^ i^^'^-^'^') (Shaver et al. 1999; Jelic et al. 2008) and the 
second comes from from the receivers. The scaling of the Aef /, the 
effective collecting area of the antennae with frequency, introduces 
also a frequency dependence. We also assume that the distribution 
of noise over the uv-plane at one frequency is Gaussian, in both 
the real and imaginary part of the visibilities. For a 24-tile LOFAR 
core station, the SEFD will be around 2000 Jy (77 ~ 0.5), depend- 
ing on the final design (de Bruyn et. al 2007). This means that we 



Here we use the Kronecker product identity AXB = C 
(B^ (g) A) vec (X) = vec(C). We can then solve for X, when 
(B^ ® A) is invertible (Golub & Loan 1989). This way the measurent 
equation can be transformed into a linear equation. 
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can reach a sensitivity of 520 mK at 150 MHz with 1 MHz band- 
width in one night of 4 hours of observation. Accumulating data 
from a hundred nights of observations brings this number down to 
~52 mK. We will assume a constant noise estimate of this mag- 
nitude for each pixel in the remainder of the paper. The noise is 
uncorrelated between different telescopes and different signals, so 
that we can write the noise matrix in the form: 



where n. 



N = 



{nT 







(17) 



We are now in a situation that we can start to simulate realistic 
data-sets for LOFAR. 



5 REALISTIC LOFAR-EOR-KSP DATA SIMULATIONS 

In this section, we will describe the LOFAR instrumental response 
simulations and the inversion method used to enhance the result of 
the primary calibration. The simulation consists of two steps: (i) 
the forward step, where we use a realistic model of the LOFAR re- 
sponse function, including all instrumental effects and noise to gen- 
erate simulated LOFAR EoR data, and (ii) an inversion step where 
we use the data model with realistic solutions and error range for 
the data model parameters, to obtain the underlying visibilities. In 
this subsection we shall consider the forward step. The data model 
is based on the Hamaker-Bregman-Sault measurement equation, 
as described in the previous sections, and we use an "onion-layer" 
approach to predict the instrumental response, as is suggested by 
the order of the Jones matrices in the Measurement Equation. First, 
we must obtain a representation of the Fourier Transform of the 
sky, as it is perceived by an interferometer. To achieve this we pre- 
dict the values of the K Jones matrix. This includes calculating the 
u, V and w terms as well as the I, m and n direction cosines. For 
the purpose of this paper we calculated them for six hours of inte- 
gration and ten seconds of averaging. The centre of the sky maps is 
at 5c = 52° and ac = 12 hours. 

The next step is to predict the 'true' imcorrupted visibilities for 
different sky realizations. Those are the visibilities as seen by a per- 
fectly calibrated and noiseless instrument. In this step, the source 
structure and the calibration problem are linked together in the 
visibilities (Pearson & Readhead 1984; Hamaker 2006, 2000a,b, 
1999), thus we need to investigate the performance of the inversion 
method, using sky models with different complexities (see subsec- 
tion 5.3). In this paper we apply the method to the galactic diffuse 
emission superimposed on the cosmological signal. In the future we 
will also consider a collection of point sources (this will be useful 
for the construction and improvement the Local Sky Model). 

The simulated maps between 200 and 110 MHz were cre- 
ated at a high frequency resolution (100 kHz) corresponding to the 
LOFAR EoR data-set resolution and were subsequently biimed in 
frequency using a 1 MHz moving average box function. Each sim- 
ulated visibility data-set has a size of approximately 100 MB. We 
produced data-sets for 128 frequency channels (reduced from 320 
that will be used during the real experiment). 

The simulation, inversion and signal extraction algorithms are 
implemented in MATLAB, with extensive use of MATLAB exe- 
cutablcs files, to speed up the most often visited loops. Due to the 
parallel nature of the procedures, the Distributing Computing Tool- 
box is used, while many BLAS operations and the EFT are imple- 
mented using the relevant NVIDIA CUDA 2.0 libraries, and run on 



a single Tesla S870 system. For the rest we used a 16-way SMP 
workstation with 32 GB of RAM and sufficient disk space. 

5.1 Cosmological signal and foregrounds simulations 

We start by discussing the cosmological signal maps that we used 
in the simulation. The cosmological 21 -cm signal is generated in 
a simplified manner, but sufficient for our purposes. We generate 
a contiguous cube of dark matter density from redshift z ~ 6 to 
« ~ 12 (corresponding to frequencies of around 110 MHz to 200 
MHz) as described in Thomas et al. (2008), from equally spaced 
outputs in time of an N-body simulation. We employ a comoving 
100 h^^ Mpc , 256'' particles dark matter only simulation, using 
GADGET 2. We furthermore assume a flat ACDM universe and 
set the cosmological parameters to f2m = 0.238, f2b = 0.0418, 
Q.K = 0.762, (78 = 0.74, = 0.951 and h = 0.73, in agreement 
with the WMAP3 observations (Spergel et al. 2007). 

The ratio between the baryons (basically atomic hydrogen and 
helium) to dark matter was set to fib/ dm ~ 0.2. Atomic hydro- 
gen is assumed to follow an analytical ionization history of expo- 
nential nature as 1 -|- exp{z — Zion)^^, where Zion differs for each 
line-of-sight as 10 ± A'^(0, 1), where A''(0, 1) is a normal distri- 
bution with mean zero and dispersion one. This is done to mimic 
a possible spread in redshift during which the reionization occurs. 
We plan to use a more realistic cosmological signal in the future, 
but for our current purposes (i.e. testing the inversion process), this 
is sufficiently complex. 

For each of the foreground components, a 5° x 5° map is 
generated in the same frequency range (between 110 and 200 MHz) 
pertaining to the LOFAR-EoR experiment. In this paper, we use 
only simulations of the Galactic diffuse synchrotron emission and 
of radio galaxies. More complex foreground realizations are under 
investigation and we will be considering them in future work. 

Of all components of the foregrotmds, galactic diffuse syn- 
chrotron emission (GDSE) is by far the most dominant and orig- 
inates from the interaction between free electrons in the interstel- 
lar medium and the Galactic magnetic field. The intensity of the 
synchrotron emission can be expressed in terms of the brightness 
temperature Tb and its spectrum is close to a featureless power law 
T(, ~ i/*^, where /3 is the brightness temperature spectral index. 

At high Galactic latitudes the minimum brightness tempera- 
ture of the GDSE is about 20 K at 325 MHz with variations of the 
order of 2 per cent on scales from 5-30 arcmin across the sky (de 
Bruyn et. al 1998). At the same Galactic latitudes, the temperature 
spectral index (3 of the GDSE is about -2.55 at 100 MHz (Shaver 
et al. 1999; Rogers & Bowman 2008) and steepens towards higher 
frequencies, but also gradually changes with position on the sky 
(e.g. Reich & Reich 1988; Platania et al. 1998). 

In the four dimensional (three spatial and one frequency) 
simulations of the GDSE, all its observed characteristics are in- 
cluded: spatial and frequency variations of brightness temperature 
and spectral index, as well as brightaess temperature variations 
along the line of sight. The spatial distribution of the 3D ampli- 
tude and brightness temperature spectral index of the GDSE are 
generated as Gaussian random fields, while along the frequency 
direction we assume a power law. The final map of the GDSE at 
each frequency is obtained by integrating the 4D cube along one 
spatial direction. All spatial and frequency properties of the sim- 
ulated GDSE are normalized to match the values of observations. 
In addition to the simulations of the total brightness temperature, 
polarized Galactic synchrotron emission maps are also produced, 
that include multiple Faraday screens along the line of sight. For a 



The WEAR Data Model 



11 



detailed explanation of the GDSE emission simulations in total and 
polarized brightness temperature see Jelic et al. (2008). 



5.2 The uncorrupted visibilities 

The prediction of the visibilities for simple intensity distribution 
like point-sources and Gaussian distributions is straightforward. 
This is not the case for diffuse emission with complex structure. 
There are two approaches in this case. In the discrete version of 
the problem, one can either use a very high resolution map and 
then assume that caeh pixel is a point source, taking into account 
the linearity of the ME over the brightness distribution, or one can 
find a proper representation of the sky distribution (i.e. shapelets) 
and use a carefully chosen gridding-convolution function on the uv- 
plane. Image plane effects need to be computed locally, as they af- 
fect relatively small scales and thus, they have a small convolution 
footprint. This is the approach of the MeqTree/UVBrick software 
(Smimov & Noordam 2006) which is used for predicting uv-data 
of extended sources or of multiple point sources. Despite the com- 
putational gains of this method, we shall use the more direct and 
precise multiple point-source prediction. This is done in order to 
preclude the introduction of spurious distortions on the simulated 
cosmological signal as well as to retain the high-resolution of the 
image-plane effect predictions. The sky brightness distribution can 
be split into patches in order to increase the speed of execution 
at the expense of additional phase-shifts (i.e. a shift in (/, m, n), 
the phase center of the image, correponds to a phase shift of the 
Fourier kernel). In that case each discrete Fourier Transform and 
each frequency channel is independent, which yields an embarrass- 
ingly parallel computational process. 

In our matrix formulation the operation performed becomes 

iVp.x 

Vtrue.f iu,V,w) = ^ KijBijKlj, 
i 

where Bi j is coherency the matrix of pixel i on a grid at each fre- 
quency /, and Kj,/ is the corresponding Fourier kernel in the di- 
rection of pixel i. The same approach holds for the point-sources as 
well, which can have any given position. Hence, we can incorporate 
the different components of our sky model into a single formalism. 
Another reason for our preference for this method, is that the inver- 
sion of the above equation is equivalent to a deconvolution process. 
We shall discuss this in more detail in the following subsection. 



5.3 Simulations of the instrumental response 

In this section, we describe the different effects the instrument and 
ionosphere have on the observed visibilities. These can be split in 
three categories: image-plane effects, uv-plane effects and noise. 



5.3.1 Image-plane effects 

As we saw, the array response matrix A is composed of multiple 
effects, namely: 

• Faraday rotation 

• Ionospheric phase fluctuations 

• Antenna voltage pattern 

• Polarization leakage and instrumental polarization 

• Parallactic angle rotation 

• Fourier kernel 



The effects that depend on s = (/, m, n) are commonly called 
image-plane effects. We apply them directly on the images be- 
fore taking into account the K Jones matrix. This is because they 
are direction-dependent gains or rotations of the coherency matrix. 
They transform each pixel on the maps in a different way. In the 
uv-plane, they enter as convolutions of the visibilities. For now, 
we ignore the parallactic angle P Jones matrix, because its terms 
should be known to very high precision a priori, and depend only 
on the rotation of the sky inside the station primary beam pattern. 

We only consider the beam pattern and the ionospheric Fara- 
day rotation and phase delays. We use the LOFAR HBA dipole 
beam patterns, as well as the HBA station patterns provided in an 
ASTRON^ technical report, by Yatawatta (2007). They depend on 
the element geometry and the station layout. Initial engineering re- 
sults from Brentjens, Yatawatta and Wijnholds (unpublished) in- 
dicate that they are stable over a time interval of ~8 hours. We 
normaUze the response at the phase center to one. Slight misalign- 
ment of the dipoles, being not precisely orthogonal, would lead to a 
different polarized beam pattern response. Because we are not cor- 
relating each dipole pair independently, but we are actually doing 
aperture synthesis between entire stations, we do not have the in- 
formation needed to correct this effect on a per dipole basis. In the 
simulations we assume that this effect is included in the errors of 
the relevant Jones matrices. We therefore include the net effect in 
the errors of the complex beam pattern. 

The other two image-plane effects that we shall consider in 
this paper are the ionospheric phase and the Faraday rotation. We 
assume that the ionosphere consists of a single layer at an altitude 
of 250 km above the ground. We assume a velocity with a mean of 
200 km and a variance of 10 km h~^. This is done in order 
to simulate the turbulent temporary motion of the traveUng iono- 
spheric disturbances (TID). We use five directions with different 
velocities. In such case the vertical and the slant Total Electron 
Contents (TEC) are equal. We assume that the TEC distribution is a 
random fields consistent with Kolmogorov-type turbulence. It has a 
mean value of 20 TECUs (1 TECU = lO" electrons per m^) and a 
variance of ~1 TECU as in Fig 5. Each station sees approximately 
20x20 square kilometers of the ionosphere, depending on the fre- 
quency. We simulated a larger patch, though, of 50 x 50 square kilo- 
meters. The ionospheric phase-delay introduced would be 



Arion 



40.3m3s-2 • TEa, 



(18) 



where TECsrc is the TEC along the line of sight, c is the speed of 
light and f is the observing frequency (Spoelstra & Yang 1990). 
The form of the ionospheric delay Jones matrix is then 

Z = Ie''^^'°° 



and the Ionospheric Faraday rotation is 



RM 



sin 



(19) 



where RM is the rotation measure (Brentjens & de Bruyn 2005). 
The rotation measure in SI units is defined as 



RM = 



http://www.astron.nl 



Jo 



ne-B|| ds 



(20) 
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where eo is the vacuum permittivity, with _B| | the magnetic field, 
m the mass of the electron, e the charge of the electron and rie 
the electron density. From the simulated maps we construct the 
coherency matrix for every pixel of the maps. We then apply the 
image-plane effects on the maps and apply the K Jones matrix to 
bring them to the uv-plane. 



5.3.2 uv-plane ejfects 

The next category of instrumental effects is the uv-plane effects: 

• Complex gains. 

• Frequency bandpasses. 

The LOFAR gains and especially the bandpasses are expected to be 
stable over the period of one night with temporal variations of the 
order of two per cent. The bandpass response is well known and 
also well behaved. For this reason, the bandpass and complex gain 
effects together are handled in a single Jones matrix G. Actually, 
the bandpasses are the frequency dependent parts of the instrumen- 
tal complex gains. The frequency dependent gain can be approxi- 
mated as: 



G(/,i) = 









G2 



G{i,2} =fl'{i,2} (l + lO ^sin(wctsi)) [l+^{v-vo)], 

(21) 

where S'{i,2} are the complex gain coefficients, uicts is the cyclic 
frequency that corresponds to the correlation time scale of the gain 
solutions, 7 is constant with a small value and vo = 150 MHz. 

Since there are no HBA station cross-correlations available at 
the moment, we used the WSRT radio telescope solutions, scaled 
up as a function of collecting area. For that we used the gain solu- 
tions from a deep survey of a few fields carried out by the LOFAR- 
EoR KSP team during November 2007 (Bemardi et al., in prepa- 
ration). The goals of this WSRT survey are to better understand 
the Galactic foregrounds, at frequencies in a range similar to that 
spanned by the LOFAR high-band antennas, and to test the calibra- 
tion pipeline. The survey was carried during 6x 12h runs. One of 
the fields observed was the so-called Fan field (Bemardi, in prepa- 
ration). Finally, we model the polarization leakage following Equa- 
tion 13, with d = and ip = O.ln this paper we ignore the 
radio-frequency interference (RFI) contamination and we assume 
that the data have been sufficiently corrected. We will consider the 
implications of insufficient RFI flagging and its effects on the sig- 
nal extraction in the future. 

Our final measurement equation, including all effects de- 
scribed above, can thus be written as: 

A{v, t) = G {v, t) D {u) E {u, t) K {u, t) F {u, t) Z {u, t) (22) 

with the parameters being defined as described in this section. 



5.3.3 Additive Noise 

Understanding the noise characteristics is very important for the 
LOFAR EoR-KSP, which aims to detect very weak radio emis- 
sion. The noise properties of the instrument are vital for a sensible 
error analysis. The total effective collecting area for the LOFAR- 
EoR experiment is ~ 0.07 km^ at 150 MHz. The instantaneous 
bandwidth of the LOFAR telescope is 32 MHz and the aim for 
the LOFAR-EoR experiment is to observe in the frequency range 



between 115 and 180 MHz, which is twice the instantaneous band- 
width. To overcome this, multiplexing in time has to be used (for 
more details see, de Bruyn et. al 2007). For the purpose of this sim- 
ulation we ignore this complication and assume 400 hours of inte- 
gration time on a single field for a bandwidth of 32 MHz, centred at 
uo = 150 MHz. This is chosen for two reasons: first, the frequency 
of 150 MHzis the frequency where reionization presumably peaks 
in the current cosmological simulations. Second, due to the size of 
the data generated and the constraints of our current hardware (a 
16-way symmetric multiprocessor system with 32 GB of memory), 
we are required to reduce the number of channels. This deteriorates 
the foreground fitting efficiency, but this is a reasonable compro- 
mise to be made. As we will show in the next section, this issue 
does not seem to pose any significant threat to our ability to statis- 
tically detect the EoR signal and our approach can thus be regarded 
as conservative. 

The ultimate sensitivity of a receiving system is determined 
principally by the system noise. The noise properties of an elab- 
orate receiving system like LOFAR can be very complex and we 
will address this more in forthcoming papers. The theoretical rms 
noise level in terms of the real and imaginary part of the complex 
visibilities for an interferometer pair between stations p and q is 
given by : 



AF{Re,Im}^, = lx 

'Is 



SEFDp X SEFDq 

2 X Al/ X Tavg 



where rjs is the system efficiency that accounts for electronic, dig- 
ital losses, Au is the frequency bandwidth and Tavg is the aver- 
aging time during which each station accumulates data. SEED is 
the System Equivalent Flux Density. The SEED can be written as 
SEFD = Tsys/K, where K = {-q^ x A^s)/{2 x fee) and de- 
pends on the station efficiency rja, and the effective collecting area 
^eff of the station, fee is the Boltzmann constant. For the system 
noise we assume two contributions. The first comes from the sky 
and is frequency dependent (« v'"^'^^) and the second comes from 
the receivers. For the LOFAR HBA stations (24 tiles) the SEFD is 
~ 2000 Jy at 150 MHz, depending on the final design (de Bruyn 
et. al 2007). We assume that the SEFD varies within one per cent 
between different stations. In order to calculate the SEFD we use 
the following system temperature (Tsys) scaling relation as function 
of frequency {v): 

Tsy, = 140 + 60(1^/300 MHz)-^-^^ 

For our simulation, we scale the noise, so that in one night of ob- 
servations we have the same sensitivity as the LOFAR-EoR ob- 
servations after 400 hours. This is done because we currently lack 
the computational capacity to simulate and analyse a full data-cube 
of order several TBs per field. The equation above gives the noise 
on the real and imaginary parts of the visibility. When the SNR 
is high the noise distribution of the visibility amplitude and phase 
is Gaussian to an extremely good approximation. Without any sig- 
nal the amplitude noise follows the Rayleigh distribution. When 
the SNR is low the measured amplitude noise follows the Rice dis- 
tribution (Thompson, Moran, & Swenson 2001; Taylor, Carilli & 
Perley 1999; Lopez & Fabregas 2002) and gives a biased estimate 
of the true underlying signal. For the phases, the noise distribution 
becomes uniform. 

In the following section we describe one method that demon- 
strates our ability to statistically detect the EoR signal from the 
post-calibration LOFAR maps that include realistic levels for the 
noise and zero-mean calibration residuals. In both cases for the sta- 
tistical detection of the signal we use the total intensity maps only. 
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Figure 5. Left: A simulated map of the TEC of the ionosphere at a particulat time, above the LOFAR core. Each LOFAR station sees approximately twenty 
square kilometers of the ionosphere at an altitude of 250 km. The values are displayed in TECUs. Right: Example of phase fluctuations caused by the 
ionosphere at a baseline of 100 metres as a function of Hour Angle 




Figure 6. A plot of the real part of the complex gain for an XX correlation 
of a 500 m baseline at 150 MHz as function of time. The dots represent the 
recovered solution while the solid line the values used to corrupt the data. 
The relative error is plotted in the bottom panel. 



In contrast to previous work, however, because we have taken po- 
larization into account, its effects leaking into Q, U and V maps 
and are therefore accounted in our analysis. Gain calibration is as- 
sumed to be done with a precision of two per cent, which is a real- 
istic estimate, judging from the experience gained thus far from the 
LOFAR Core Station One (CSl) data analysis. 



5.4 The data-model inversion method 

To invert the ME, we have rewritten it as a linear equation, relating 
the observed visibilities to the true underlying visibilities, using the 
properties of the outer Kronecker product (Golub & Loan 1989; 
Boonstra 2005). The relevant Jones-Mueller matrix describes the 
parameters of the model. This way we have to deal with a linear 
(over the parameters) model. Each parameter is a non-linear, multi- 
variate function of other parameters that describe the physics be- 
hind each of the effects described by the relevant Jones matrices. 



Due to our lack of information, though, we shall consider the ele- 
ments of the Jones-Mueller matrix as multilinear functions of those 
parameters. Our task is to infer the parameters of the Jones-Mueller 
matrix, that we sample in the presence of Gaussian calibration er- 
rors and white noise. 

The maximum-likelihood (ML) is a powerful method for find- 
ing the free parameters of the model to provide a good fit. The 
method was pioneered by R. A. Fisher (Fisher 1922) . The maxi- 
mum likelihood estimator (MLE) selects the parameter value which 
gives the observed data the largest possible probability density in 
the absence of a prior, although the latter can be easily incorporated. 
For small numbers of samples, the bias of maximum likelihood es- 
timators can be substantial, but for fairly weak regularity conditions 
it can be considered asymptotically optimal (Mackay 2003). With 
large numbers of data points, such as in the case of the LOFAR 
EoR KSP, the bias of the method tends to zero. In general, it is not 
feasible to estimate the size of the data needed in order to obtain a 
good enough degree of approximation of the likelihood function to 
a multivariate Gaussian. We assume that errors associated with each 
datum are independent, but they have different variances and corre- 
lation properties over time and/or frequency. Many of these prob- 
lems can be overcome by MCMC (Markov Chain Monte Carlo) or 
nested sampling of the posterior (Skilling 2004; Feroz, Hobson, & 
Bridges 2008). 

Traditionally, the CLEAN/self-calibration approach has been 
used in radio-interferometry to give estimates of both the model pa- 
rameters and the sky brightness distribution (Pearson & Readhead 
1984). Given an initial model of the sky distribution, which is usu- 
ally made using CLEAN components for a few bright sources with 
a known structure, we solve for the non-linear parameters p of the 
model, and then improve the initial map. This is done in a loop until 
the optimization converges. In our approach, we make a distinction 
between the calibration and the data-model inversion parts. This is 
justified by the volume of the data to be handled. Another factor 
which is relevant for the EoR KSP is that we have to work near 
and below the noise level, and initial estimates of the sky distribu- 
tion are going to be affected by the noise and calibration errors. We 
therefore need to make a proper determination of the true underly- 
ing visibilities and the model parameters, while retaining the noise 
and cosmological signal properties. 

In our approach we assume that the parameter vector p is 
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known within an accuracy of one to two per cent. For example, 
in Fig. 6, the gain used in the simulation and the recovered solution 
are shown. The parameter vector p includes all the (non)-linear pa- 
rameters of the instrumental and ionospheric models, that are part 
of the Jones matrices (see section 4.2). This is a key assumption and 
its impact needs to be determined. We will address this issue in the 
future, when the LOFAR roU-out will provide us with more data to 
test our current error estimation levels. Given that, the maximum 
likelihood solution can be formulated in a more powerful matrix 
formalism, as: 

[m] (p) N-4,Mi (p)] b. = [mJ (p) N"!,] vob.,., 

where Nnoise is the determinable noise covariances and M = 
A (g) A. Its inverse plays the role of a metric in the N-dimensional 
vector space of the data. The metric is useful for answering ques- 
tions having to do with the geometry of a vector space. In our case 
we use it to actually compute the dot product on that space. Solving 
this equation gives the ML solution for Gaussian noise. By lineariz- 
ing the above ML equation along the eigenvectors of p, we expect 
to find a better solution. Eventhough this can be done fully analyt- 
ically, the complexity of the resulting matrices is so large that we 
plan to determine the later through forthcoming numerical simu- 
lations. By perturbing p, we can then recover the proper noise and 
EoR covariance matrices and determine our final solution of s. This 
process could even be iterated, but we expect convergence after a 
couple of iterations, since we start with an already good approxi- 
mation of the model parameters. With the addition of extra regu- 
larization terms, a deconvolved image can be obtained, from which 
the foreground and point-sources can be subtracted or filtered, to 
leave only the EoR signal and the noise. This 'ideal' process, will 
be the subject of the second paper in this series. In the current pa- 
per, however, our main goal is to show that if calibration errors can 
be corrected without bias (i.e. errors on the correction have asymp- 
totically zero mean), that we can indeed recover the EoR signal. 
Clearly, this is the first step in this process. In the following section 
we will therefore apply the more simple method used by Jelic et al. 
(2008), which does not do a ML inversion, but simply uses "dirty" 
images with identical synthesized beam shapes, in order to extract 
the cosmological signal from the dirty maps. 



6 RESULTS FROM THE SIGNAL EXTRACTION 

The ultimate benchmark whether we can, in principle, recover the 
redshifted 21 -cm signal from the LOFAR EoR KSP data set, is to 
apply our ML inversion and our signal-recovery procedure to the 
simulated data with realistic sky, ionospheric and instrument set- 
tings, and if that is successful, ultimately to real data sets. 

This, however, is not yet computationally feasible^" and we 
opt in this paper for a less computationally expensive approach by 
removing the foregrounds from the "dirty" images. We are in fact 
implementing the full ML method on a mini-cluster of three quad- 
core PCs connected to three NVIDIA Graphics Processor Units. 

To remove the foregrounds, we currently use a polynomial fit 
along frequency for every Une of sight on the "dirty" images. To 
produce the dirty maps we use only uv points that are present in 
95% of the frequency channels (Figure 7). This leads to an identical 
PSF for every map, with the additional cost of resolution loss. The 

Current tests indicate 30 hrs of computing time for the ML inversion of 
a single frequency channel of realistic data size on a single high-end CPU 



Table Lvalues of the noise icTnoise [niK] at 150 MHz), for different levels 
of calibration errors 



case (a) 


case (b) 


case (c) 


case (d) 


Error [%] 0.5 


1 


2 


10 


o-noise[inK] 36 


58 


60 


73 



visibilities are gridded on a super-sampled, regular grid of 256 by 
256 cells. This degrades the resolution of the dirty maps somewhat 
but minimizes any mixture between spatial and frequency fluctua- 
tions because of gaps in the uv-plane. Our current analysis is also 
done including only the Galactic foregrounds, and we will address 
the issue of the extragalactic foregrounds in the future, once we 
have more computational power available to us. 

An inappropriate polynomial fitting, however, could remove 
part of the EoR signal or in the case of under fitting of the fore- 
grounds, fitting residuals could dominate over the EoR signal (for 
details and discussion see Jelic et al. (2008)). Hence, after sub- 
stracting the foregrounds from the data cubes, the residuals should 
ideally contain only the noise plus the EoR signal. The noise level 
is nearly an order of magnitude larger than the EoR signal, how- 
ever, so one is able to make only a statistical detection of the signal 
over the map by taking the difference between the variance of the 
residuals and the variance of the noise. The underlying assumption 
here is that the general statistical properties of the noise are known 
to a high accuracy as a function of frequency. We are currently in- 
vestigating methods that provide accurate noise estimates from the 
or partly calibrated data. 

Fig. 8 shows the standard deviation of residuals as a func- 
tion of frequency (dashed line), after taking out the smooth compo- 
nent of the foregrounds using a third-order polynomial. The white 
dashed line represents the mean of the detected EoR signal after 
100 independent Monte Carlo simulations of the extraction method 
applied to each realization, while the grey shaded zone shows 2cr 
detection limits. As one can see the detected EoR signal is in good 
agreement with the original (solid red line) for most the frequen- 
cies (Table 1). We use bold-face letters to mark our current expec- 
tations. We thus conclude that in this case we can still recover the 
EoR signal, while for larger errors we gradually lose our ability to 
statistically detect the EoR signal. 



7 SUMMARY & OUTLOOK 

The main goal of this paper has been to introduce a physics-based 
data model for the LOFAR EoR Key Science Project, use this data 
model to generate reaUstic data sets, including Galactic and ex- 
tragalactic foregrounds and known instrumental, ionospheric and 
noise properties (with zero-mean residuals), and subsequently ex- 
tract the redshifted 21 -cm EoR signal from the data cubes. De- 
spite simplifications, which we indicated and will further address 
in forthcoming papers, we have clearly made a large step forward 
in simulating realistic circumstances under which currently planned 
EoR experiments like LOFAR will operate. This has thus far been 
lacking in many publications on this topic, which have all neglected 
calibration and uv-sampling issues, and clearly needs to be ad- 
dressed beyond overly simplified assumptions before results from 
any of the currently planned projects can be believed. 

We have also shown that using the currently known or esti- 
mated instrument specifications and layout of LOFAR, the 21 -cm 
EoR signal can be recovered if the noise properties are known a 
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Figure 7. The sampling of the uv-plane by the LOFAR core along frequency after 6 hours of synthesis. The left-hand figures show the average number of 
visibilities per uv cell for 8, 16, 32 and 64 MHz of total bandwith (the instantaneous bandwith of LOFAR is 32MHz). We assume that the data are delivered at 
0. 1 MHz resolution. The colorbar shows the number of visibilities per grid point. The right-hand set of figures shows the area in which less than 5% of the data 
along frequency is lost due to the scaling of the uv coverage with frequency, compared to the the total bandwidth. The black points represent regions where the 
visibilities and their Fourier conjugates occupy the same place, while the grey points represent true visibility measurements. This distinction is made because 
the Fourier conjugates do not contribute to the SNR. 



priori (although we have not yet tested whether the extraction of 
noise and EoR signal separately is possible), and the calibration of 
the ionosphere and instrument can be carried out in an unbiased 
fashion to the level that WSRT and LOFAR test setups indicate is 
possible. In forthcoming papers we will address far more complex 
situations that include all possible foregrounds (not only Galactic) 
and where calibration will be done using the simulated data them- 
selves. This requires a considerable increase in our computational 
power and we are currently implementing our codes on GPUs to 
achieve the required processing power 

In more detail, in this paper we initially started by providing a 
brief description of the physical connections behind the Hamaker- 
Bregman-Sault formalism. This completes the picture established 
by Hamaker & Bregman (1996). The connection to physics helps 
in better understanding new and old problems in radio-polarimetry 
and provides a new perspective. Novel methods of calibration and 
image evaluation can be developed based on those principles and 
we indicated a few possibilities (see Appendix A). 

After describing the connection between the physics of the 
Jones and Mueller calculi and the data model for a generic radio 
interferometric array, we introduced the data model that is relevant 
for the LOFAR EoR KSR This description of the signal pathways 
was used to generate detailed, large-scale polarized instrumental 
response simulations for several instrumental parameters assuming 
simplified, but still realistic models for their values and errors. The 
high-resolution data produced by this simulation were generated in 
such a way that they resemble a typical LOFAR EoR observation. 

We applied the same procedure, albeit simple, as described 
by Jelic et al. (2008) to recover the cosmological signal, and as a 
first benchmark of the inversion method. The polynomial fit method 
(Zaldarriaga, Furlanetto, & Hemquist 2004; Jelic et al. 2008) gives 
a reasonable result, even in the presence of small errors in the cal- 
ibration parameters. When including realistic instrumental noise 
and calibration errors, we can still recover the cosmological signal 
if the properties of instrumental noise are well known as a func- 
tion of frequency. In order to obtain the frequency dependence of 
the noise, we used the true underlying sky distribution. In reality. 



one might not have this luxury. However, differencing of image or 
visibility values between narrow frequency channels could in prin- 
ciple deliver the noise properties of the instrument as function of 
frequency. 

What has not been addressed in this paper, but will be in the 
follow-up paper, is the calibration process itself. The traditional 
self-calibration algorithm (Pearson & Readhead 1984) alternates 
between calibration and CLEAN steps in order to find the optimal 
image, from a given visibility data set. This is done by substracting 
simple models for the sky brighness distribution and correcting for 
the calibration solutions at every step. This, however, leads to a sub- 
optimal solution and certainly cannot deal easily with complex sky 
brightness structures (Starck, Pantin, & Murtagh 2002; Bhatnagar 
& Cornwell 2004; Puetter & Pina 1994; Narayan & Nityananda 
1986; Cornwell & Evans 1985; Kemball & Martinsek 2005). On 
the other hand, the maximum likelihood method has traditionally 
been used in many disciplines in order to estimate the parameters 
of linear systems. It can address the calibration, (after lineariza- 
tion), and deconvolution problems and in forthcoming papers we 
will further test this methodology, besides others [e.g. Expectation 
Maximization; Yatawatta et al. (2008)]. 

Finally, there is the problem of the signal extraction itself. 
Although done independently of the inversion process in this pa- 
per, we suspect it can possibly be done simultaneously with the 
inversion. This will also be addressed in forthcoming papers. Here, 
we have studied signal extraction through polynomial fitting, but 
matched filtering is another approach that we plan to take, and 
we suspect it might be a good complementary method. Those al- 
gorithms might have to make some use of a priori assumptions, 
though, since the observational evidence for the cosmological sig- 
nal is scarce to non-existent. 

A detailed analysis of the noise properties of the foreground- 
cleaned maps also still has to be carried out. This can be done in 
a Monte Carlo sense. This way one can study the asymptotic effi- 
ciency of the ML, using the observed information matrix. Another 
consideration should be the correlation properties between the in- 
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Figure 8. The signal extraction for (a) 0.5 per cent (top left), (b) 1 percent (top right) , (c) 2 per cent (bottom left) and (d) 10 per cent (bottom right) errors 
on the calibration parameters. The solid blue line represents the rms value of the residuals of the fiting procedure. The black-dotted line is the rms of the 
instrumental noise as a function of frequency. The rms of the "dirty" map of the cosmological signal is plotted with the red line and of the extracted signal with 
the white dashed line. The shaded grey areas are the errors on the rms of the extracted signal calculated from 100 realizations of the noise. The effectiveness 
of the extraction proceedure decreases rapidly with the amplitude of the errors. 



strumental parameters, as well as the error propagation, since in 
principle their behaviour is highly non-linear. 

Despite a considerable increase in the complexity of large- 
scale EoR simulations (Thomas et al. 2008), foregrounds models 
(Jelic et al. 2008) and the instrument model (this paper) - necessary 
to process and understand the results from our forthcoming LO- 
FAR EoR KSP observations - with this paper we intend to provide 
a guide towards even more complex simulations, including all the 
effects that we mentioned above and thus far have avoided or im- 
plemented in a simplified manner. We conclude with the comment, 
learned from the cuiTent paper, that only by doing simulations and 
successful signal-recovery on a scale and complexity-level com- 
parable to the real observations, can one convincingly show that 
ongoing experiments - be it LOFAR, MWA, or otherwise - are 
capable of extracting the redshifted 21 -cm signal from data sets af- 
fected by noise, RFI, calibration errors and bright foregrounds. It is 
also essential in interpreting their ultimate results. 
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APPENDIX A: PROPERTIES OF THE STOKES VECTOR 
AND MATRIX 

In this appendix we introduce several fundamental concepts in the 
description of polarized radiation. Although not yet used in this 
paper, we feel that they are important to mention and will be used 
in forthcoming publications to guide us to a better understanding of 
complex caUbration processes in radio interferometry. 



Al Pauli Matrices 

Pauli matrices are well known and have been used for the analy- 
sis of partially polarized light (Fano 1954). Their major advantage 
is that they satisfy a set of properties that significantly reduce the 
complexity of calculations associated with the intensity. The iden- 
tity plus the Pauli matrices in two dimensions are defined as 



CO = 



IT2 = 



1 

1 

1 

1 



<T3 



1 

-1 

-i 

■I 



(Al) 



This set of 2 x 2 linearly independent matrices constitute a basis 
for the vector space of 2 x 2 Hermitian matrices over the com- 
pex numbers. Surrmiarizing their properties, they are Hermitian and 
they follow the commutation relations \ai, aj\ = aiaj — ajai = 
i2sijk(Tk where eijk is the Levi-Civita permutation symbol (Golub 



& Loan 1989; Fano 1954; Boonstra 2005). These matrices are uni- 
tary and traceless except for the identity matrix. The Unear expan- 
sion of the coherency matrix in this basis is 



(A2) 



with .s, = tr (C(T,;) being the four Stokes parameters: i = 0, 1, 2, 3 
corresponding to the Stokes /, Q, U and V, respectively. 

A2 The Stokes Vector and Matrix 

In the Uterature the Stokes parameters are usually arranged as a 
4x1 vector, s = {I, Q, U, V)"^ (Gil 2007). An alternative notation 
would be to introduce a 2 x 2 Stokes matrix. 



I + Q U-iV 
U + iV I-Q 



(A3) 



In this case for pure states we have ||C|1'^ = | and in the case 
of unpolarized light ||C||'^ = | ||S||^. This is because for unpolar- 
ized light Q, U and V are zero and the matrix becomes I multiplied 
by the identity matrix. Since there is a factor of |, the Frobenius 
norm of the matrix has this factor squared. Even if the source is not 
polarized one should observe the same signal in both orthogonal 
polarizations, as in this case both polarization states are equiproba- 
ble. This is an important fact for calibration. The X and Y dipoles 
of the LOFAR HBAs should measure the same signal, for unpo- 
larized sources, and any fluctuations should be due to polarization 
calibration errors. 



A3 Polarization Level 

In addition, the degree of polarization is defined as P s 
y[/2 + v'^ + Q^/I. We can also define the unit vector (Gil 2007) 




The Stokes vector can be decomposed using those parameters in 
several ways. A trivial decomposition is between a polarized and 
unpolarized state. As radio receivers are intrinsically polarized, and 
in the case of LOFAR orthogonal, a spectral decomposition can ex- 
press the Stokes vector as a convex linear sum of two orthogonal 
pure states. This makes the relationship between the measured sig- 
nals and the parameters describing the system more clear, since we 
have to deal with orthogonal states. The spectral decomposition of 
the coherency matrix is using its eigenvalue structure to decompose 
it into pure states, i.e. eigenvectors. The spectral decomposition (di- 
agonalization) of the coherency matrix is then equivalent to 



s = / X 



Of course, for a mixed state there are infinite combinations of 
independent states into which it can be decomposed, but it would 
be useful if the selection is such that it matches the characteristics 
of the system. The parameters / and P are invariant under unitary 
transformations (changes of coordinate systems). They are directly 
related to the eigenvalues of C. Pure states are related to rank-1 po- 
larization matrices, and mixed states to rank-2. Wolf showed that 
there always exist two orthogonal reference directions, such that 
the degree of coherence is maximized and coincides with P. This is 
important because various random distributions correspond to un- 
polarized light. The measurement of the correlations of the Stokes 
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parameters allows us to distinguish between those different types 
of non-polarized Ught. 

A4 Entropy 

Concluding with the quality criteria, the von Neumann entropy 
can be applied to electromagnetic waves (Fano 1957; Brosseau 
1998). In terms of the coherency matrix it is defined as S = 
— tr In . Is a measure of the difference in the amount of in- 
formation between pure and mixed states, both with same inten- 
sity. It can be expressed as a function of the eigenvalues of C as 
S = -i [(l + P)ln(yrTP) + (l-P)ln(Vr^)].Itisa 
decreasing monotonic, bounded function of P (Gil 2007). It at- 
tains its maximum value, S — ln(2), for non polarized light and 
its minimum, S — 0, for totally polarized light. Using this entropy 
one can define the polarization temperature. Depolarizing effects, 
such as the ionospheric Faraday rotation, can be studied this way 
in order to detect spatial heterogeneity. For example, ionospheric 
TIDs should lead to an increase in the polarization entropy. This is 
a good scheme of ranking observations: if the polarization entropy 
differs between two maps, residual Faraday rotation and leakage 
can be present. Finally, for non-Gaussian distributions of polariza- 
tion states, higher order moments are needed and the Stokes system 
is no longer adequate to describe the polarization states. 

The standard (self)-cahbration procedure for interferometric 
observation of polarized Ught aims at recovering the coherency ma- 
trix of the polarized radiation. Given the statistical nature of the co- 
herency matrix, we must emphasize the importance of parameters 
that give a measurement of their polarimetric purity. Polarization 
entropy is a concept related to the impurities of the media through 
which radiation propagates, as it was defined above. It is useful for 
many purposes, especially when depolarization is a relevant sub- 
ject. An increase in the polarization entropy signifies a decrease 
in the polarization purity. It is a direct way to access the quaUty 
of the data. Ionospheric Faraday rotation causes a depolarization 
with a certain frequency behavior. We would expect that the en- 
tropy would follow this behavior. In Figure Al we show the polar- 
ization entropy for several lines of sight as a function of frequency 
in a simulated map. Attention should be brought to the fact that is 
really hard to distinguish between Faraday rotation and other depo- 
larization effects. 

A5 Lorentz lyansformations 

The transformations of the Stokes vector lead to a system of differ- 
ential equations (Landi Degl'Innocenti 1987). From the statistical 
interpretation of the coherency matrix we derive that the Stokes I 
parameter has to be positive and that so = I ^ sf + S2 + sf^ — 
+ + V^ which means that there can be no more polarized 
light than the total light. This is a key concept in physics, namely 
the conservation of energy. Extending this remark, we can define 
a Minkowski space (in the mathematical sense) for the Stokes 4- 
vectors. The norm in this space would be ||s|| = I^ — Q^+U^ + V^. 
Positive values of the norm correspond to the time-like vectors of 
the special theory of relativity. Light-like Stokes vectors correspond 
to totally polarized states. The Poincare sphere defines the light 
cone with the exception that the symmetric part of the cone does 
not have any special meaning. The Jones vector transforms as usual 
under the Lorentz group. The coherency matrix is defined through 
the Kronecker product of two Jones vectors and is a 2x2 spinor 
of rank two. The Stokes vector and the coherency matrix must be 



realizations of the same irreducible representation of the Lorentz 
group, as they both have 4 independent components. The Lorentz 
transformations form a 6 parameter group. The six generators are 
the 3 spatial rotation generators of 0(3) R and the Lorentz boosts 
B. The infinitesimal transformation of the Stokes vector is 
3 

s' = s — (nRi + 6iBi)sd7 = s — Ksd7 

i 

where K is a matrix that resembles the absorption matrix. This 
equation is the radiative transfer equation along the signal path. If 
one choses to exclude effects that intensify the radiation, this is the 
complete mathematical description of the effects in the signal path. 
The Lorentz boosts describe polarizing effects, while the spatial 
rotations describe the Faraday rotation. The relation between the 
initial and final Stokes vector is a finite Lorentz transformation. 
While it sounds straightforward it is not an easy task. The problem 
arises from the non commutativity of the generators. Magnus has 
proposed a solution to this problem. We will discuss this soon. 

A6 Clifford Algebra 

After having introduced the coherency matrix and the Jones for- 
malism we proceed a step further in the mathematical abstraction 
in order to unify all possible cases. We note that every system is 
equivalent to a parallel combination of pure systems, and any pure 
system to a combination of retarders and diattenuators. In Hamaker 
(2000c) the author proposes the use of the quaternion algebra to 
describe the Jones matrices, or their 50^(1,3) covering group. 
Quaternions, as any other space isomorphic to C constitute a sub- 
algebra of a Clifford algebra. In particular, in three dimensions the 
Clifford product is equivalent to spatial rotations. Pure systems can 
still be described as Lorentz transformations in this algebra. Since 
algebraic manipulations become more clear, the computational re- 
quirements of problems involving partial polarization can be re- 
duced. The Clifford algebra of three-dimensional space CZ3 rep- 
resents the four dimensional Minkowski space time. Clifford al- 
gebras have been used extensively to describe polarization mode 
dispersion (PMD) in optical fibres (Reimer & Yevich 2006). In that 
field they have to deal with a spatially inhomogeneous, birefringent 
medium (the optical fibre), which resembles the effects caused by 
the ionosphere on the EM waves that propagate through it. 

A7 Magnus Expansion 

In the case of polarization mode dispersion, which can occur due 
to ionospheric Faraday rotation, transversal of the polarized source 
signal through a set of phase screens etc, one is interested in re- 
covering the original polarization state of the radiation. Two possi- 
bly quasi-orthogonal modes, represented by the input frequency- 
independent and the output Jones vectors, are related through a 
complex 2x2 Jones matrix. Despite the fact that we might not have 
any prior knowledge about the structure of the medium through 
which radiation propagates, the frequency dependence of that ef- 
fect contains useful information about that medium. For example 
the ~ Faraday rotation should leave a distinct imprint on the 
signal, which can help us distinguish it from other depolarizing ef- 
fects. In the general case, the coherency matrix is transformed as 
C = JCJ^ An arbitrary Jones matrix with determinant 1, can be 
also expressed in terms of two vectors a and b according to 

J = exp [- (i/2) (b + ia) ■ a] , (A4) 




where a is the Pauh spin vector. This is a Jones matrix representa- 
tion of the Lorentz transformation in a Clifford algebra. In the case 
of a frequency dependent effect, 3 {uj), the Jones space operator 
can be decomposed into real and imaginary components as 

J {lo) = exp [- (j/2) (b + ia) • a] , 

- - i (A5) 
J^J = --[n + iA], 

where the subscript cu represents differentiation with respect to the 
frequency. A solution for J can be obtained through the Magnus 

expanion, which specifies J (uj) = exp I ^ Bn (ui) J J (luo), 

\ri=0 / 

where the first two coefficients are given by 

(^) = £^ dujiduji (J^ (wi) J„ (wa) - {(^2) (oJi)). 

The coefficients of the Magnus expansion for n > 2, are re- 
lated to those of lower order through recursion. Taylor expanding 
the frequency derivatives of the Jones matrix to third order gives 
us a way to directly evaluate the Magnus coefficients. To the extent 
of our knowledge, the methods mentioned above have not been ap- 
plied in the field of radio interferometry. 



